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Abstract 

A combination of implicit and explicit timestepping is analyzed for 
a system of ODEs motivated by ones arising from spatial discretiza- 
tions of evolutionary partial differential equations. Loosely speaking, 
the method we consider is implicit in local and stabilizing terms in 
the underlying PDE and explicit in nonlocal and unstabilizing terms. 
Unconditional stability and convergence of the numerical scheme are 



proven by the energy method and by algebraic techniques. This sta- 
biUty result is surprising because usually when different methods are 
combined, the stability properties of the least stable method plays a 
determining role in the combination. 

1 Introduction 

This report considers timestepping methods for systems of ordinary differen- 
tial equations of the form 

u'{t) + Au{t) + B{u)u{t) - Cu{t) = f{t), (1.1) 

in which A, B{u), and C are nxn matrices, u{t) and f{t) are n- vectors, and 

A = A^yO, B{u) = -B{u)\ C = ^ and A - C ^ 0. (1.2) 

Here >- and ^ denote, respectively, the positive definite and the positive 
semidefinite ordering. The key properties motivating our work are that A is 
sparse and that although C is not sparse, the action of C on a vector is inex- 
pensive to calculate. This structure is motivated by multiscale discretizations 
of turbulence but can also arise from closed-loop control problems and en- 
semble calculations. Given this structure of (jl.ip . the simplest scheme that 
is computationally feasible is explicit in the global, unstable part of ()1.1|1 . 
that is, Cu. Accordingly, we consider 

— + AUn+l + B{Un)Un+l - Cu^ = fn+1, k = At, (1.3) 

k 

where Un is the approximation to u{t = nk). Usually when methods are 
combined, the stability properties of the explicit method play a determining 
role in the overall method. In Theorems 12 . II and 12 . 21 we prove the surprising 
result that ()1.3|1 is unconditionally stable. This result is outside the realm of 
root condition stability analysis for uncoupled scalar problems. 

In Section 2, unconditional stability and convergence of ()1.3|1 are proven. 
We give two stability proofs. The first is algebraic. Since the constants 
depend on the dimension of the system, we also give an energy proof of 
stability (with uniform constants) that is potentially extensible to discretized 
PDEs. Section 3 presents numerical tests illustrating the theory. First, we 
briefly summarize some motivating problems leading to (jl.lll . 
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The basic model of the turbulent dispersion is that it is dissipative in the 
mean (see [TUj, [Ej, [7j). A more accurate formulation is that its dissipative 
effects are focused on the smallest resolved scales (see |5]). This physical 
idea has led to algorithms for numerical stabilization of transport-dominated 
phenomena based on eddy diffusivity acting only on the smallest resolved 
scales (e.g., 0, 0, HSI, 0, 0, H, 0, CH, US)- The natural realization 
of this idea for spatial discretizations of convection diffusion equations is 
diffusive stabilization on all scales and then antidiffusing on the large scales. 
This leads to the system of ODEs 

iiijit) + b ■ V^'uij - (eoih) + e)A''uij + eo{h)PH{A''PH{uij)) = fi,, (1.4) 

where standard notation is used: A'* is the discrete Laplacian, eo{h) is the 
artificial viscosity parameters and Ph denotes a projection onto a coarser 
mesh; see Section 3 for details. The system (jl.4p fits exactly the form (jl.ip . 
(II. 2j) . where C is provided as the matrix arising from eo(/i) term. We shall 
also test one algorithm as a perturbation of the method ()1.4|) in which the 
projection is replaced by a nearest averaging A^UiJ. In both cases, the projec- 
tion or averaging operator accounts for the nonlocal character (i.e., the large 
bandwidth) of C. On the other hand, averaging and projection are both 
embarrassingly parallel operators whose action on a given vector is cheap to 
perform. 



Remark 1.1. (1) A second main application is discretization of turbulent 

flow problems which, although nonlinear and constrained, have a similar 

structure to the above (simple) linear convection diffusion problem. 

(2) A known method of stabilizing the timestepping and the associated linear 

system (but not the spacial discretization ) corresponds to / li.i)) without the 

averaging: 

!h^±l-^ + b ■ V"n„+i - (eoih) + e) A^„+i + eo(/i)A^„ = (1.5) 

Each time step requires the inversion of the matrix corresponding the operator 
— (eo(/i) + e)A'^ + 6- V'^ + Zc"^/, which, fore^ suitably chosen, is an M-matrix. 
Our analysis applies to this method as well. 
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2 The Stability Analysis 



For our analysis, we assume that B{u) is in C^ijif^) and /(t) is in C"'^([0, oo)). 
For any T > 0, we denote by 



pt = niax II / 1 1. Ml 9 . 

t€[0,T] 

Lemma 2.1. The system of ODEs M.l]) under the condition il.^) with initial 
condition u{0) = uq has a unique solution on [0, T], for any T > 0. 

Proof Since p.lj) can be written as m = i/jlt, u) with if) being of class C° 
in t and in u, local existence and uniqueness follows from the classical 
theory of ODEs j2i Theorem V.8]. 

We now show that the solution does not experience blow-up and can be 
extended everywhere. We multiply through by {u{t)'^) and we use (|1.2|) 
to obtain that 

u{tfu\t) < -u{tf{A - C)u{t) + u{tff{t) < u{tff{t). 
Using Cauchy-Schwarz, we obtain that 



j^\\u{t)\\l<\\u{t)\\l + Fl 



In turn, this implies that, 

||n(t)||^<|h(0)||^e* + F|(e*-l) 

for any t in an interval containing where u{t) is defined. Since u{t) does 
not experience blow-up in finite time, it can be extended uniquely over all of 
[0,T]. □ 

Note that from and our assumtion that f{t) is of class C^([0, cx)]), 
we get that u{t) is of class C^([0, oo]). The fact that u"{t) is continuous will 
be used in determining a bound for the truncation error. 

Consider the system of ODEs p.lj] under the condition p.2|l and dis- 
cretized by ()1.3p . 

First, we note that each step of ()1.3|) requires the inversion of I+kA+kBn- 

Lemma 2.2. Under / li.jjj) the n x n matrix I + kA + kB^ has a positive 
definite symmetric part and is invertible. 



4 



Proof: Let x be any nonzero vector in 3?". Then 

x^{I + kA + kBn)x = x'^x + kx^Ax + kx'^BnX 

= llxW^"^ + kx^Ax > 0. □ 

Since A, = B{un) and C do not commute, the stabihty of the numeri- 
cal scheme cannot be analyzed by reduction to eigenvalues. Therefore, we 
formulate an energy norm that is not increased at each time step, that is, 

||lfn+l lis — ll'^'"ll E- 

Definition 2.1. The energy norm of M.'J^) . is given by 

||w||g^ = u^u + ku^Cu, (2.1) 

for some u G 3?", and its associated inner product is < u,v >e= (Nkv)'^ (Nku) , 
with Nk = {I + kC)2 J for some -u, f G SfJ". 

It can be seem immediately that the energy norm and the 2-norm satisfy 
the following inequality: 

a/1 + kXmin{C) < \\u\\ ^ < + kXmax{C), 

where Xmin{C) and Xmax{C) are, respectively, the smallest and the largest 
eigenvalue of C. From this inequality and the positive semidefiniteness of C, 
we get that the induced matrix norms satisfy 

\\A\\e < \\A\\^ ^/l + kX„^ax{C). 

Theorem 2.1. Let Un satisfy MJ^) with /(.) = 0, under the condition M/<3^) 
on the coefficients. Then, 

ll'^n+lll E — ll^"ll E ■ 

Proof: Multiplying with through the equation in ()1.3p . we obtain 

Un+1 ^ r '^n+l^''^n+l + Un+l-'^nUn+l — U^^iLyUn 

Since Bn is skew symmetric, u^j^^BnUn+i = 0. Therefore 

Un+1 + U^.+lAUn+l = U^+iCUn. (2.2) 
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This is equivalent to 
Since C, we have that 

< ul^^Un + kul^^Cu^. (2.4) 

Define w = {un+i, k^^^C^^^Un+if ,v = {un^k^l^C^'^UnY ■ Then (Q) can 
be written as w'^w < w^v. Applying the Cauchy-Schwarz inequality, we get 
||w||2 < ||'i^||2- Hence, 

U^^-^Un+l + ku^j^-^CUn+l < U^Un + ku^CUn (2.5) 

or 

The conclusion of the preceding theorem is that when (jl.ip is homoge- 
neous, / = 0, we obtain that llMnlU 7^ ll^oL) ^'^j independent of T. This 
means that our method is, indeed, unconditionally stable. 

Consider (jl.3|) with / = 0, rewritten as 

{I + kA + kBn)Un+l = {! + kC)Un, S„ = B{Un) . (2.6) 

Equation (|2.6|) yields 

Un+i = {I + kA + kBn)-\l + kC)un, 

which, in turn, implies that 

(/ + A;C)5m„+i = (/ + A:C)5(/ + kA + kBnY\l + kC)"^{I + kCf^Un. 

Therefore, from the definition of |H|^, a sufficient condition to prove the 
unconditional stability result is to prove that 



(/ + kC)-^{I + kA + kBn)-\l + kC) 
From fT!^ this can be done by using the following Lemma. 



< 1. 
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Lemma 2.3. Let Di = Di"^ >- and D2 = D-i^ y be n x n matrices 

1 

such that Di — D2 >~ 0. Let D4 = D| and be symmetric. If is an n x n 
skew-symmetric matrix, then 

\\D^{Di + Ds)-'DY\2<l. (2.7) 
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Proof: Let F = D^{Di + D^) ^D^. It is straightforward that 
= D^^{Di + Di)D^^. For any nonzero vector x in S?", 

x^F-^x = x^Dl^{Di + D^)D^^x 

= x^D-^DiD^^x + x'^D-^DsD-^x 

Here we claim that D^^ D-^D^^ is skew symmetric and therefore x^ D^^ D-^D^^x = 
0. To obtain this one can notice that since D2 = D\ and D2 is a symmetric 
matrix, then D4 and D^^ are also symmetric. 
Hence, 

Thus 

x'^F-^x = x^D^^DiDl^x, for any ^ z G W. 
Using the fact that Di — D2 is nonnegative, we obtain 

x^F-^x > x^D^^D2Dl^x = x^x, for any ^ x G 3?". 

This implies that 

||2;||2^ < x^F^^x < \\x\\^ . \\F^^x\\^ , for any 7^ x G W^, 

that is, 

11^112 < ll^^^^lls' for any ^ X G 3?". (2.8) 
Obviously ()2.8|) is equivalent to 

\\Fy\\,< \\y\\,, for any ^ y G 3?" (2.9) 

Since the last equation holds for any nonzero vector y, then ||F||2 < 1. 
□ 

For the next step, we analyze the stability of the nonhomogenous problem 
over an arbitrary but finite time interval [0,r]. We later show that the 
stabihty of the homogeneous problem does not depend on T. Consider ()1.3p 
with / ^ 0. 

After some simple calculations, we get that ii„ satisfies 

Un+i = {I + kA + kBn)-\l + kC)un + k{I + kA + kBny^fn+1. (2.10) 

We denote the range of the step index n, by [0,A^], where kN = T. To 
simplify the notation, we do not explicitly indicate that depends on k and 
T. 
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Theorem 2.2. Let hold. Then the solution of \2.10\} satisfies the fol- 
lowing hound: 

k 

^ 11^°"-+ H + kLm ll^(^) II. , VO < n < ^ - 1. 
Here T is the size of the integration interval. 

Proof: To simplify notation, we take Nk = {I + kC)^ and Mk = (/ + kA + 
kBn)'^{I + kC). Then the equation ()2.1()j) can be written as 

Un+l = MkUn + k{I + kA + kBn)~^fn+l- 

Using the definition 2.1, we have 

{NkUn+lf{NkUn+l) = {NkUn+lfNkMkUn + k{NkUn+lfNk{I+kA+kBny^fn+l. 

Algebraic manipulation and the Cauchy-Schwarz inequality yield 



||A^fc^Xn+l||2 < \\NkUn+l\\^.\\NkMkNj;'\L.\\NkU 



n\\2 



+ k ||iVfcn„+i||2 • \\NkMkN^\ . ||A^fc-Vn+l||2 • 

Using Lemma EHl with D2 = and Di + D3 = MkN^'^, we obtain that 
\\NkMkN,^\\^ < 1. Then the previous inequality reduces to 



k 112 

\\NkUr,+i\\^ < ||iVfcM„||2 + k ||iVfc"Vn+l| 



This inequality can be simplified as follows: 

||A^fcM„+i|l2 < \\NkUn\\2 + k\\N-^Nkfn+il^ 

< \\NkUj, + k\\{I + kC)-'\\^\\Nkf n+l\ 



k 

(1 + kXmin{C)) 



^ ll^fc«n|l2 + 7^^^ 77^\\Nkfn+l\\2 



Thus, 

k 



Un+l\\E hn\\E< (1 + fcA^^^ (C) ) " ^"+1 " ^ 
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and since (/ + kC) ^ is a symmetric positive definite matrix, 

1 



|/ + A;C||2 = maxA(/ + A;C) 



-1 



(minA(J + A;C))' 

By the spectral mapping theorem A(/ + kC) = 1 + k\{C). Therefore 
where Xmin{C) is the minimum eigenvalue of matrix C. This implies 



Summing from to n gives 

k 

Ikn+ilU - \\uo\\e < (i + A;A ■ (C)) ^^=° < n < - 1, 

that is, 

k 

Wuu+iWe < \\uo\\e + (i + A;A ■ ((g)) ^^=° ll-^P+i"^^ ' < ^ < ^ - 1> 

which is the claimed first result. The second result follows immediately. □ 

The local truncation error of the method (jl.3|) is clearly 0{At). In the 
error estimate (which follows) we need a precise statement of this fact, which 
we now derive. 

To simplify our notation, we use Un to denote u{tn), where u{-) is the exact 
solution of (jl.ip . We also use Un to denote an iterate of our numerical scheme, 
but the particular meaning of Un will become evident from the context. 
According to the definition of local truncation error pP , 

k 

-[u'{tn+l) + Au{tn+l) + B{u{tn+l))u{tn+l) - Cu{tn+l)] (2.11) 
= ^"^^^ - u'^+l - {B{u{tn,+l)) - B{u{tn)))u{tn+l) + C(u„+i - U„). 

Using the second-order integral form of the Taylor expansion around tn+i, 
we obtain 

Un+1 - Un - ku'^j^^ = - U"{t){t - t„+i)rft, 
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which we rewrite as 



k 



- - <+i = -T / U"{t){t - tn+l)dt = -T / U"{t){tn+1 - t)dt. 
Jt„+i Jt„ 



Using the first-order integral form of the Taylor expansion around tn, we 
obtain 

- Cu'{t)) dt. 

Using the expression we have derived for the local truncation error Tn+i, 
and the preceding equations derived from Taylor's theorem, we obtain 



r, 



n+l 



~kJt 

Jtn \ 



U"{t){tn+l-t)dt- 



/ d 
di 



B{u{t))u{tn+i) - Cu'it) dt 



^-U"{t) - ^B{u{t))u{tn+l) + Cu\t) ) dt. 



k 



dt 



By the mean value theorem, there exists ^„ e (tn, tn+i) such that 



Tn+i = -U"{^n){tn+1 -Q-k J^BHt)) 



uitr,+i) + kCu'iU). (2.12) 



Hence, using the fact that < (t„+i — < k, we obtain that 



||r„+i||2 < A; max (| 

tn<S<tn+l y 



t=S 



rnax \\u{e)\\, + \\Cu'{s)\\, 



This proves the following lemma. 

Lemma 2.4. Let k — At and n>0. The method 

■^71+ 1 ~ '^n 



k 



+ AUn-\-i + B^Un+l — CUn — fn+1) 



(2.13) 



where A — >- and C — >z are n x n symmetric matrices, Bn an 
n X n skew-symmetric matrix, and fn+i — f{{'n+ l)k), is consistent. That 
is, the local truncation error is 0{At). 

We now bound the total error. We consider first the energy norm of 
truncation error. 
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Lemma 2.5. Let r„+i be the local truncation error of method \2.1'J\) . Then 

d 



I'^n+illp < k max I \\u" 



E + W^^'i'^JWE 



dt 



B{u{t)) 



max uisl p 



(2.14) 



Proof: By definition of energy norm and following the identity ()2.12|) . we 

get 



-u'\in){tn+l~in)-kj^B{u{t)) 



U{tn+l) + kCu\in) 



t=^n 



for some € [tn,tn+i]- The conclusion follows after applying the inequality 
< tn+i — Cn ^ k, the triangle inequality, and the properties of the max 
function. Note that is the induced ||-||^ of the corresponding 

matrix. □ 



We now give a convergence result for the solution of (jl.3|) . First , we need 
to compute a certain estimate. We have that 

d 

- = J —[B{u{tn)9 + Un{l - d))]u{tn+l)dt = 

{Vu{B{u{tn)9n + Mn(l " ^n)))e„) for SOmC dn G [0, 1], 



where e„ = u{tn) — Un- Here u{tn) is the solution of ()1.1|) . whereas u„ is the 
solution of our numerical scheme. 

We define the matrix VF„, by its action on a vector x G 3?": 

WnX = [{WuBiu{tn)9^ + U„(l - 9^))) x] w(t„+i), 

which results in the following identity 

- B{Un)] U{t^^,) = WnCn. (2.15) 

Lemma 2.6. Let u(.) he the solution of M-1]) and Un he the approximation 
to u{nAt), ohtained from the numerical scheme ( li.,^) . Then there exists T 
such that, Vt G [0, T] we have that 

\\Wn\\^ < r, and <rE = T^l + kXmax{C), \/0<n<N. 
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Proof: From Theorem 12.21 we have that 



u 



E 



< VI + k\max{C) (ll^olla + Tmaxj6[o,T] II/WII2) , VO < n < A^. 
We define 

Ai? = Vl + TA™ax.(C) (\\u4^ + T max ||/(t)||2 ) • 

From Lemma l2.ll we have that u{t) is bounded on [0,T], and we define 
Ku = maxtg[o,T] ||'^(0ll2- Since B{-) is of class C^, we can define 

r= max \\[{\/uB{9v2 + il-9)ui))u2]v,\\^. 

From the definition of Wn, we immediately obtain that 

||W„||2 < F, VO < n < A^. 

The second part of the conclusion follows from the inequality between H-H^, 
and ll'llg. □ 

Theorem 2.3. Consider solving the nonhomogenous problem on the interval 
[0,TJ 

u' + Au + B{u)u -Cu = f 
using the following method 



k 



+ + BnUn+l — CUn = fn+1, 



where k = At, Bn = B{un) and fn+i = f{{n + 1)A;). Let e„ = u{tn) — Un 
denote the local error. Assume that cq = 0. Then the method is convergent 
and 

)n— 1 
~ ^ k'^U 

\^n+l\\E ^ r~ Wl^ 



1 -|- kXuii'fi 



ei+fc^mmC^) — 1 k'^U 

l-\-kXYnin (C) 
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< kT 77^, VO < n < A^- 1, 



when r^; 7^ 0, and 

\\en+l\\E < ('^ + 1) 

when = 0, where 



< T- 



kU 



1 ~|- kXjyiiji (c)- 1 + kx^uc) 



\/0 < n < N -1, 



0<t<T 



U=max{\\u"{t)\\^ + \\Cu'it)\\^ + 



max mis) p 



Proof: Following the definition of the truncation error t„_|_i and using the 
equation (j2.15p . we obtain that the error, e„ = u{tn) — Un, satisfies 



k 



+ ACn+l + BnCn+l - C Cn = Tn+l - WnCn- 



After algebraic calculations, we find that 

Cn+i = {I + kA + kBn)'\l + kC)en + k{I + kA + A;S„)-^(r„+i - WnCn). 
We use the energy inner product to obtain 

< Cn+i, e-n+l >E = 

<{I + kA + kBn)-\I + kC)en + k{I + kA + A;S„)-1(t„+i - Vr„e„), e„+i >e ■ 

Applying the definition of energy norm ()2.1|) and the substitutions Mk = 
{I + kA + kBn)~\l + kC), and A^^ = (/ + fcC)5, we find that 

[NkCn+iY NkMkCn + k{Nken+iYNk{I + kA + A;S„)-1)(t„+i - iy„e„). 
Using the Cauchy-Schwarz inequality, we obtain that 



+k\\Nken+ih.\\NkMkN-^ 



2 

-1 1 



--nil 2 
r-1 



Thus 



n+l II2 



Using Lemma 12.31 with D2 = A| and Di + = M^N^ ^, we obtain that 
< 1. Hence 



|AfcMfcA,-i| 



|e„+l||^ < ||e„||^ + k llA'fc Wn+l - WnCn^E^ 
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Equivalently, we obtain that 

||e„+i||£; < ||e„||^ + A;||(/ + A;C)"^||2(||r„+i||^+ || W^„||^ ||e„||^) . 

Notice that (/ + kC)~^ is a symmetric positive definite matrix and 

1 



\{I + kC)-^ 



1 -|- kXffiin 



On the other hand, by LemmaEEl there is a constant Te such that || W^H^; < 
Te- Therefore, 

/ kVE A II II ^ 

This is a recursion formula of the following form: 

rn+i < + 6r„, 
which, when a 7^ has an upper bound of the type 

a^~^ - 1 

Tn+i < a'^ro H 6max||7:„||^ . 

a " 

Using this fact, we obtain that, when Te ^ 0, the following bound for 
the error holds whenever < n < — 1. 



1 + ^ ) _1 , 

J- ^ l+fcA™,„(C) J k ,1 ,1 

Replacing ||r„+i||^ by its bound ()2.14p obtained in Lemma and con- 
sidering that Co = 0, we have, when Te ^ and < n < — 1, that 

|en+l|| E — 
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with U = maxo<t<T {\\u"{t)\\^ + \\Cu'{t)\\^ + \\f^B{u{t))\\^m8iXo<s<T \Hs)\\e). 
The second inequahty for F 7^ follows from the inequality (1 + a;)" < e^", 
for X > and n positive integer. 

When r^; = 0, we immediately get from ()2.16|) and from Lemma f2. 51 that 

which, together with kN = T prove the inequalities for Te = 0. 

The convergence follows from the fact that ||-||^ converges to ||-||2 as 
A; — >■ which implies that ||e„||2 — as /c — 0. □ 

The case F^; = occurs, for example, when B{u) is constant (which we 
simulate numerically in the next section). For that case, the error increases 
only linearily with the size of the interval, assuming that the derivatives up 
to order 2 of the solution u{t) are uniformly bounded. 



3 Numerical Results 

Let Q = [0, 1] X [0, 1]. For the equation 

Ut + b ■ Vu — eAu 
u 

u{x, 0) 



= /, over 

= (j){x) on 6n, (3.1) 
= ^0(2;) in Q, 



use the method described in this work, with uniform mesh and central dif- 
ference. A choice must be made for the antidiffusion operator: averaging 
or projection. We have selected averaging. Since it is just outside the the- 
ory, we will thereby test the robustness of the algorithm. Antidiffusion is 
completed by averaging, where u (]?):= weighted average of nearest neighbors. 
This corresponds to filtering with 6 = 2h. The method becomes in our case 

Uij{t) + b ■ V^Uij — (e + eo)A^Mjj + eoAuij'^ = h, 

where q denotes how many times the average operation is taken. In our 
experiments we chose q = 2 and e = 10~^. We take b = (cos(6'), sin(6')), 
where 6 = 17°. 

For the boundary and initial conditions we take the line at angle 9 through 
the center of the domain. On the north side of the line we take (p = 1 on the 
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£g=0.0001 



£p=0.001 




25 



Energy norm for ^^^=0. 0001 




I time step = 0.1 

time step = 1 

,1 1 \ \ \ \ \ \ I I I 

01 23456789 10 

time 



Figure 2: Stability of the numerical method demonstrated by the behavior 
of the energy norm 

Energy norm and shifterj energy norm for e ^=0.005 




0123456789 10 

time 



Figure 3: Numerical validation of Theorem 12.11 

17 



Explicit advection treatment, timestep=0. 1 , e ^=0.005 



10 - 




01 23456789 10 

time 

Figure 4: Exponential growth of the solution of the scheme that includes 
the advection term explicitly 

boundary; on the south side we take = on the boundary. We take / = 
and as initial conditions. 

We performed the following experiments, all on a 32 x 32 mesh. 

1. We ran the simulation for 1,000 steps with a timestep of 10, with 
the artificial viscosity parameter eo having succesively the values 10~^, 
5 X 10^^, 10^^, 10"'^. We have presented no analysis for the spatial 
dependence of the solution with respect to eo, but we have included 
this experiment for validation, since our choice of parameters should 
result roughly in the steady-state approximation for this mesh, which 
has been studied before in the literature. 

The results are depicted in Figure ^ We see that when the artificial 
viscosity parameter eo is very small, a complete loss of coherence of the 
spatial structure results, whereas too large a parameter (cq = 0.1) alters 
the steady-state solution significantly. This effect is consistent with the 
typical behavior of centered methods for the skew step problem |T3] . 

2. For cq = 10"'^, we ran the simulation for 100 steps with a timestep of 1 
and for 1, 000 steps with a timestep of 0.1. The energy norm comparison 
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of these computations is presented in Figure |21 We see that even for 
the very large step, the energy norm stays bounded, consistent with 
our absolute stability claim. 

We also present in Figure El a comparison between the energy norms of 
the distance between the successive iterates of the two cases and their 
outcome at time 100. From Figure |21 we infer that m(100) is a rea- 
sonable approximation to the steady-state solution. Since the equation 
i\'S.l\i is linear, we have that Un — m(100) is the result of the numeri- 
cal scheme applied to the homogeneous equation associated to ()3.1|) . 
From Theorem 12. II we have that — it(100)||£; must be a decreasing 
sequence, which is exactly what we observe from Figure El Note that 
ll-u^ll^ is not a decreasing sequence, as can be seen in Figure El More- 
over, the sequence H^nll^; may not even be monotonic, as seen in Figure 
H for A: = 0.1. 

3. We compare the results of our scheme with the similar scheme that 
takes into account explicitly the term that contains the skew-symmetric 
matrix B{un). For the latter scheme we obtain the recursion 

— + AUn+l + B{Un)Un - CUn = fn+1- 

K 

We apply this scheme to our example on a 32 x 32 mesh for 1000 
timesteps of length k = 1. We see the rapid exponential growth that 
is typical for computations with the timestep outside the region of 
stability. 

This demonstrates that our scheme has significantly better stability 
properties than the alternative, which would result in linear systems of 
comparable sparsity. The numerical scheme, based on a backward Euler 
approach that considers all terms implicitly, though absolutely stable, 
will result in less sparse linear systems since the matrix C contains 
an averaging operator that substantially reduces sparsity and is not 
considered here for comparison. 
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